#include "PSFRun.hh"
#include "G4AnalysisManager.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void PSFRun::Merge(const G4Run *run)
{
  const auto localRun = static_cast<const PSFRun *>(run);

  fNTIR += localRun->fNTIR;
  fNTIR2 += localRun->fNTIR2;
  fNExiting += localRun->fNExiting;
  fNExiting2 += localRun->fNExiting2;
  fEscapedEnd += localRun->fEscapedEnd;
  fEscapedEnd2 += localRun->fEscapedEnd2;
  fEscapedMid += localRun->fEscapedMid;
  fEscapedMid2 += localRun->fEscapedMid2;
  fBounce += localRun->fBounce;
  fBounce2 += localRun->fBounce2;
  fPSFCoreBounce += localRun->fPSFCoreBounce;
  fPSFCoreBounce2 += localRun->fPSFCoreBounce2;
  fClad1Bounce += localRun->fClad1Bounce;
  fClad1Bounce2 += localRun->fClad1Bounce2;
  fClad2Bounce += localRun->fClad2Bounce;
  fClad2Bounce2 += localRun->fClad2Bounce2;
  fReflected += localRun->fReflected;
  fReflected2 += localRun->fReflected2;
  fEscaped += localRun->fEscaped;
  fEscaped2 += localRun->fEscaped2;
  fMirror += localRun->fMirror;
  fMirror2 += localRun->fMirror2;
  fDetectorHits += localRun->fDetectorHits;
  fDetectorHits2 += localRun->fDetectorHits2;
  fCoreTIR += localRun->fCoreTIR;
  fCoreTIR2 += localRun->fCoreTIR2;
  fCoreR += localRun->fCoreR;
  fCoreR2 += localRun->fCoreR2;
  fClad1TIR += localRun->fClad1TIR;
  fClad1TIR2 += localRun->fClad1TIR2;
  fClad1R += localRun->fClad1R;
  fClad1R2 += localRun->fClad1R2;
  fCoupleR += localRun->fCoupleR;
  fCoupleR2 += localRun->fCoupleR2;

  G4Run::Merge(run);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void PSFRun::EndOfRun()
{

  if (numberOfEvent == 0)
    return;
  auto TotNbofEvents = G4double(numberOfEvent);

  fNTIR = fNTIR / TotNbofEvents;
  fNTIR2 = fNTIR2 / TotNbofEvents;
  G4double rmsTIR = fNTIR2 - fNTIR * fNTIR;
  if (rmsTIR > 0.)
    rmsTIR = std::sqrt(rmsTIR);
  else
    rmsTIR = 0.;

  fNExiting = fNExiting / TotNbofEvents;
  fNExiting2 = fNExiting2 / TotNbofEvents;
  G4double rmsExiting = fNExiting2 - fNExiting * fNExiting;
  if (rmsExiting > 0.)
    rmsExiting = std::sqrt(rmsExiting);
  else
    rmsExiting = 0.;

  fEscapedEnd = fEscapedEnd / TotNbofEvents;
  fEscapedEnd2 = fEscapedEnd2 / TotNbofEvents;
  G4double rmsEscapedEnd = fEscapedEnd2 - fEscapedEnd * fEscapedEnd;
  if (rmsEscapedEnd > 0.)
    rmsEscapedEnd = std::sqrt(rmsEscapedEnd);
  else
    rmsEscapedEnd = 0.;

  fEscapedMid = fEscapedMid / TotNbofEvents;
  fEscapedMid2 = fEscapedMid2 / TotNbofEvents;
  G4double rmsEscapedMid = fEscapedMid2 - fEscapedMid * fEscapedMid;
  if (rmsEscapedMid > 0.)
    rmsEscapedMid = std::sqrt(rmsEscapedMid);
  else
    rmsEscapedMid = 0.;

  fBounce = fBounce / TotNbofEvents;
  fBounce2 = fBounce2 / TotNbofEvents;
  G4double rmsBounce = fBounce2 - fBounce * fBounce;
  if (rmsBounce > 0.)
    rmsBounce = std::sqrt(rmsBounce);
  else
    rmsBounce = 0.;

  fPSFCoreBounce = fPSFCoreBounce / TotNbofEvents;
  fPSFCoreBounce2 = fPSFCoreBounce2 / TotNbofEvents;
  G4double rmsPSFCoreBounce = fPSFCoreBounce2 - fPSFCoreBounce * fPSFCoreBounce;
  if (rmsPSFCoreBounce > 0.)
    rmsPSFCoreBounce = std::sqrt(rmsPSFCoreBounce);
  else
    rmsPSFCoreBounce = 0.;

  fClad1Bounce = fClad1Bounce / TotNbofEvents;
  fClad1Bounce2 = fClad1Bounce2 / TotNbofEvents;
  G4double rmsClad1Bounce = fClad1Bounce2 - fClad1Bounce * fClad1Bounce;
  if (rmsClad1Bounce > 0.)
    rmsClad1Bounce = std::sqrt(rmsClad1Bounce);
  else
    rmsClad1Bounce = 0.;

  fClad2Bounce = fClad2Bounce / TotNbofEvents;
  fClad2Bounce2 = fClad2Bounce2 / TotNbofEvents;
  G4double rmsClad2Bounce = fClad2Bounce2 - fClad2Bounce * fClad2Bounce;
  if (rmsClad2Bounce > 0.)
    rmsClad2Bounce = std::sqrt(rmsClad2Bounce);
  else
    rmsClad2Bounce = 0.;

  fReflected = fReflected / TotNbofEvents;
  fReflected2 = fReflected2 / TotNbofEvents;
  G4double rmsReflected = fReflected2 - fReflected * fReflected;
  if (rmsReflected > 0.)
    rmsReflected = std::sqrt(rmsReflected);
  else
    rmsReflected = 0.;

  fEscaped = fEscaped / TotNbofEvents;
  fEscaped2 = fEscaped2 / TotNbofEvents;
  G4double rmsEscaped = fEscaped2 - fEscaped * fEscaped;
  if (rmsEscaped > 0.)
    rmsEscaped = std::sqrt(rmsEscaped);
  else
    rmsEscaped = 0.;

  fMirror = fMirror / TotNbofEvents;
  fMirror2 = fMirror2 / TotNbofEvents;
  G4double rmsMirror = fMirror2 - fMirror * fMirror;
  if (rmsMirror > 0.)
    rmsMirror = std::sqrt(rmsMirror);
  else
    rmsMirror = 0.;

  fDetectorHits = fDetectorHits / TotNbofEvents;
  fDetectorHits2 = fDetectorHits2 / TotNbofEvents;
  G4double rmsDetectorHits = fDetectorHits2 - fDetectorHits * fDetectorHits;
  if (rmsDetectorHits > 0.)
    rmsDetectorHits = std::sqrt(rmsDetectorHits);
  else
    rmsDetectorHits = 0.;

  fCoreTIR = fCoreTIR / TotNbofEvents;
  fCoreTIR2 = fCoreTIR2 / TotNbofEvents;
  G4double rmsCoreTIR = fCoreTIR2 - fCoreTIR * fCoreTIR;
  if (rmsCoreTIR > 0.)
    rmsCoreTIR = std::sqrt(rmsCoreTIR);
  else
    rmsCoreTIR = 0.;

  fCoreR = fCoreR / TotNbofEvents;
  fCoreR2 = fCoreR2 / TotNbofEvents;
  G4double rmsCoreR = fCoreR2 - fCoreR * fCoreR;
  if (rmsCoreR > 0.)
    rmsCoreR = std::sqrt(rmsCoreR);
  else
    rmsCoreR = 0.;

  fClad1TIR = fClad1TIR / TotNbofEvents;
  fClad1TIR2 = fClad1TIR2 / TotNbofEvents;
  G4double rmsClad1TIR = fClad1TIR2 - fClad1TIR * fClad1TIR;
  if (rmsClad1TIR > 0.)
    rmsClad1TIR = std::sqrt(rmsClad1TIR);
  else
    rmsClad1TIR = 0.;

  fClad1R = fClad1R / TotNbofEvents;
  fClad1R2 = fClad1R2 / TotNbofEvents;
  G4double rmsClad1R = fClad1R2 - fClad1R * fClad1R;
  if (rmsClad1R > 0.)
    rmsClad1R = std::sqrt(rmsClad1R);
  else
    rmsClad1R = 0.;

  fCoupleR = fCoupleR / TotNbofEvents;
  fCoupleR2 = fCoupleR2 / TotNbofEvents;
  G4double rmsCoupleR = fCoupleR2 - fCoupleR * fCoupleR;
  if (rmsCoupleR > 0.)
    rmsCoupleR = std::sqrt(rmsCoupleR);
  else
    rmsCoupleR = 0.;
  //////////////////////////////////////////////

  G4int prec = G4cout.precision(3);
  G4cout << "\n ======================== run summary =====================\n";
  G4cout << "Average number per event:" << G4endl;
  G4cout << " TIR:          " << fNTIR << " +- " << rmsTIR << G4endl
         << " Exiting:      " << fNExiting << " +- " << rmsExiting << G4endl
         << " Escaped Mid:  " << fEscapedMid << " +- " << rmsEscapedMid
         << G4endl << " Escaped End:  " << fEscapedEnd << " +- "
         << rmsEscapedEnd << G4endl << " Bounced:      " << fBounce << " +- "
         << rmsBounce << G4endl << " PSFCore Bounce:" << fPSFCoreBounce << " +- "
         << rmsPSFCoreBounce << G4endl << " Clad1 Bounce: " << fClad1Bounce
         << " +- " << rmsClad1Bounce << G4endl
         << " Clad2 Bounce: " << fClad2Bounce << " +- " << rmsClad2Bounce
         << G4endl << " Reflected:    " << fReflected << " +- " << rmsReflected
         << G4endl << " Escaped:      " << fEscaped << " +- " << rmsEscaped
         << G4endl << " Mirror:       " << fMirror << " +- " << rmsMirror
         << G4endl << " Detector hit: " << fDetectorHits << " +- " << rmsDetectorHits
         << G4endl << " CoreTIR:      " << fCoreTIR << " +- " << rmsCoreTIR
         << G4endl << " CoreR:        " << fCoreR << " +- " << rmsCoreR
         << G4endl << " Clad1TIR:     " << fClad1TIR << " +- " << rmsClad1TIR
         << G4endl << " Clad1R:       " << fClad1R << " +- " << rmsClad1R
         << G4endl << " CoupleR:      " << fCoupleR << " +- " << rmsCoupleR
         << G4endl;
  G4cout << G4endl;
  G4cout.precision(prec);
}